# STAT 415/615 Regression (M. Baron)

# Python Lab 5. REGRESSION DIAGNOSTICS

# Run regression and look at residual plots. Split the plot window to get all the plots at once.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy import stats

# STUDENTIZED RESIDUALS AND OUTLIERS

# Load the mtcars data.
# statsmodels includes this data set.

mtcars = sm.datasets.get_rdataset("mtcars").data

# Fit the regression model.

X = sm.add_constant(mtcars["wt"])
y = mtcars["mpg"]
reg = sm.OLS(y, X).fit()

# Plot the data and regression line.

plt.scatter(mtcars["wt"], mtcars["mpg"])
plt.plot(mtcars["wt"], reg.predict(X))
plt.xlabel("wt")
plt.ylabel("mpg")
plt.show()

# Look at the standard regression diagnostic plots.

fig = plt.figure(figsize=(10, 8))
sm.graphics.plot_regress_exog(reg, "wt", fig=fig)
plt.show()

# Studentized residuals and testing for outliers

# Studentized residuals

influence = reg.get_influence()
t = influence.resid_studentized_external

# Return to the 1x1 plot window

plt.figure()

# See if there are any nonlinear trends

plt.plot(t, "o")
plt.xlabel("Observation")
plt.ylabel("Studentized residual")
plt.axhline(0)
plt.show()

# Find residuals whose absolute value is greater than 2.

t[np.abs(t) > 2]

array([2.32816206, 2.53780106, 2.38384376])

# Which of these residuals can be considered as outliers?

# Compare with the Bonferroni-adjusted quantile from t-distribution.

n = len(mtcars["wt"])

qt = stats.t.ppf(0.025 / n, n - 2)
qt

-3.478735536774798

# Find residuals exceeding the Bonferroni-adjusted cutoff.

t[np.abs(t) > np.abs(qt)]

array([], dtype=float64)

# The test revealed no outliers.

# Testing NORMALITY

# Shapiro-Wilk test

stats.shapiro(t)

ShapiroResult(statistic=0.929157838596398, pvalue=0.037105119804681015)

# Rather marginal. Also look at the Normal Q-Q plot among residual plots above.

# It is not straight, so the data may be non-Normal. Shapiro-Wilk statistic W
# measures how close the graph is to a straight line.

# Testing HOMOSCEDASTICITY (constant variance).

# Look at the residual plots including a plot of t^2 vs fitted values.

fitted = reg.fittedvalues
plt.scatter(fitted, t**2)
plt.xlabel("Fitted values")
plt.ylabel("Studentized residuals squared")
plt.show()

# Below is the Breusch-Pagan test for non-constant variance.

# The test is available in statsmodels.

from statsmodels.stats.diagnostic import het_breuschpagan

bp_test = het_breuschpagan(reg.resid, reg.model.exog)

# The output contains several statistics.

# The first two are the Breusch-Pagan statistic and its p-value.

bp_test

(0.04043839551874129,
 0.8406258536969659,
 0.03795896453698885,
 0.8468390237860595)

# This package also has a built-in outlier test that calculates the
# Bonferroni-adjusted p-values for the externally studentized residuals.

studentized = influence.resid_studentized_external
p_values = 2 * stats.t.sf(np.abs(studentized), reg.df_resid)

bonferroni_p = np.minimum(p_values * n, 1)

# Display the observation with the largest absolute studentized residual.

largest = np.argmax(np.abs(studentized))

studentized[largest], p_values[largest], bonferroni_p[largest]

(2.5378010594812555, 0.01658652213251693, 0.5307687082405418)

# No Studentized residuals with Bonferonni p < 0.05

# LACK OF FIT TEST

# The ToothGrowth dataset has only 3 different values of X = dose

ToothGrowth = sm.datasets.get_rdataset("ToothGrowth").data

ToothGrowth.head()

# Check the variable names.

ToothGrowth.columns

Index(['len', 'supp', 'dose'], dtype='object')

# Count the observations at each value of dose.

ToothGrowth["dose"].value_counts().sort_index()

# Fit two regression models:

# reduced = simple linear regression predicting Y = length in terms of X = dose

X = sm.add_constant(ToothGrowth["dose"])
reduced = sm.OLS(ToothGrowth["len"], X).fit()

# full = using group means to predict Y for each value of X,
# thus treating X as a categorical variable

dose_dummies = pd.get_dummies(ToothGrowth["dose"], drop_first=True, dtype=float)
X_full = sm.add_constant(dose_dummies)
full = sm.OLS(ToothGrowth["len"], X_full).fit()

# Plot the data.

plt.scatter(ToothGrowth["dose"], ToothGrowth["len"])

# Plot the reduced model in red.

plt.plot(ToothGrowth["dose"], reduced.predict(X), linewidth=4)

# Plot the predictions from the full model in blue.

plt.scatter(ToothGrowth["dose"], full.predict(X_full), linewidth=10)

plt.xlabel("dose")
plt.ylabel("len")
plt.show()

# Here is the rigorous F-test for the lack of fit.

# Compare the residual sums of squares from the two models.

RSS_reduced = np.sum(reduced.resid**2)
RSS_full = np.sum(full.resid**2)

df_reduced = reduced.df_resid
df_full = full.df_resid

df_difference = df_reduced - df_full

F = ((RSS_reduced - RSS_full) / df_difference) / (RSS_full / df_full)

p_value = stats.f.sf(F, df_difference, df_full)

F, p_value

(11.231909566634299, 0.0014321769516278155)

# Conclusion: the difference in SSreg is significant; the linear regression model
# does have a lack of fit.

# Box-Cox Transformation

# Find the best power transformation of responses that fixes non-normality.
# Fit a linear regression model, save studentized residuals, and test their
# Normal distribution.

X = sm.add_constant(mtcars["wt"])
reg = sm.OLS(mtcars["mpg"], X).fit()

influence = reg.get_influence()
t = influence.resid_studentized_external

stats.shapiro(t)

ShapiroResult(statistic=0.929157838596398, pvalue=0.037105119804681015)

# Results are marginal, so let’s look for the best transformation.

# scipy.stats has a Box-Cox function.

# Box-Cox requires positive response values.

# We search for the value of lambda that maximizes the likelihood.

lambdas = np.arange(-2, 2.01, 0.01)

log_likelihood = []

for lam in lambdas:
if lam == 0:
Z = np.log(mtcars["mpg"])
else:
Z = (mtcars["mpg"]**lam - 1) / lam

```
model = sm.OLS(Z, X).fit()
log_likelihood.append(-len(Z) / 2 * np.log(np.sum(model.resid**2) / len(Z)))
```

plt.plot(lambdas, log_likelihood)
plt.xlabel("lambda")
plt.ylabel("Log likelihood")
plt.show()

# Following the maximum likelihood principle, we are looking for the value of

# lambda that maximizes the likelihood function, or equivalently, the logarithm

# of this likelihood, which is on the graph.

# We see that the optimal lambda is somewhere between -1 and 0.

# We can limit the search to this range.

lambdas = np.arange(-1, 0.01, 0.01)

log_likelihood = []

for lam in lambdas:
    Z = (mtcars["mpg"]**lam - 1) / lam
    model = sm.OLS(Z, X).fit()
    log_likelihood.append(-len(Z) / 2 * np.log(np.sum(model.resid**2) / len(Z)))

plt.plot(lambdas, log_likelihood)
plt.xlabel("lambda")
plt.ylabel("Log likelihood")
plt.show()

# Now we can see that the best lambda is very close to -0.2.

# Let’s introduce a variable that is the corresponding power transform of our
# response Y, fit this new regression, and check residuals for Normality.

Z = mtcars["mpg"]**(-2)

newreg = sm.OLS(Z, X).fit()

influence = newreg.get_influence()
t = influence.resid_studentized_external

stats.shapiro(t)

ShapiroResult(statistic=0.972585398441951, pvalue=0.5736137412259258)

Z = mtcars["mpg"]**(-0.2)

newreg = sm.OLS(Z, X).fit()

influence = newreg.get_influence()
t = influence.resid_studentized_external

stats.shapiro(t)

ShapiroResult(statistic=0.9589239733012132, pvalue=0.2566072020434485)

# This regression model passes the Shapiro-Wilk test for Normality;

# the p-value is high. No evidence that this assumption is violated.

	len	supp	dose
0	4.2	VC	0.5
1	11.5	VC	0.5
2	7.3	VC	0.5
3	5.8	VC	0.5
4	6.4	VC	0.5